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Abstract 

We present the results of global three dimensional magneto-hydrodynamic simu- 
lations of black hole accretion flows. We focus on the dependence of numerical results 
on the gas temperature Tout supplied from the outer region. General relativistic effects 
are taken into account using the pseudo-Newtonian potential. We ignore the radia- 
tive cooling of the accreting gas. The initial state is a torus whose density maximum 
is at 35rs or 50rs from the gravitating center, where Vs is the Schwarzschild radius. 
The torus is initially threaded by a weak azimuthal magnetic field. We found that 
mass accretion rate and the mass outflow rate strongly depend on the temperature of 
the initial torus. The ratio of the average Maxwell stress generated by the magneto- 
rotational instability (MRI) to gas pressure, a = {B^B^/in) / (P) is a ~ 0.05 in the 
hot torus (Tout ~ 9.5 x lO^K at 50rs) and a ~ 0.01 in the cool torus (Tout ~ 1.1 x lO^K 
at 35rs). In the cool model, a constant angular momentum inner torus is formed 
around 4 — Sr^. This inner torus deforms itself from a circle to a crescent quasi- 
periodically. During this deformation, the mass accretion rate, the magnetic energy 
and the Maxwell stress increase. As the magnetic energy is released, the inner torus 
returns to a circular shape and starts the next cycle. 

Power spectral density (PSD) of the time variation of the mass accretion rate in 
the cool model has a low frequency peak around lOHz when we assumed a IOMq 
black hole. The PSD of the hot model is fiat in 1 — 30Hz. The slope of the PSD in 
the cool model is steeper than that in the hot model in 30 — lOOHz. The mass outflow 
rate in the low temperature model also shows quasi-periodic oscillation. Intermittent 
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outflows are created in both models. The outflow speed is 0.01c — 0.05c. The mass 
outflow rate is comparable to the mass accretion rate M at 2.5rs in the high temper- 
ature model and about 0.02M in the low temperature model. 

Key words: accretion, accretion disks — black hole physics — 
magnetohydrodynamicsiMHD — QPO — dynamo 

1. Introduction 

Accretion disks are believed to be the energy source of various activities in active galactic 
nuclei, X-ray binaries, protostars and so on. Using the RXTE satellite, detailed timing analysis 
of X-ray binaries has been carried out (e.g., Homan et al. 2005, McClintock & Remillard 
2006, Shaposhnikov et al. 2006). Black hole candidates (BHCs) show transitions between 
a low/hard state (LHS) dominated by the hard power-law component and a high/soft state 
(HSS) characterized by the soft black body component. During the transition, BHCs stay 
in hard intermediate state (HIMS) or in a soft intermediate state (SIMS) (Homan & Belloni 
2005). The light curves during LHS of BHCs are subject to violent X-ray fluctuations and 
sometimes accompany quasi-periodic oscillations (QPO) in the Fourier Power Spectral Density 
(PSD). The LHS is associated with steady outflows (e.g., Gallo et al. 2003). In this state, 
mass accretes to the black hole as an optically thin, advection-dominated accretion flow (e.g., 
Ichimaru 1977; Narayan & Yi 1994, 1995). The energy spectrum of the HIMS is softer than 
that in the LHS. Low frequency (1 — lOHz) QPOs are observed in the luminous hard state and 
in the intermediate state (McClintock & Remillard 2006). High frequency QPOs (~ lOOHz) are 
sometimes observed in these states. The frequency of low-frequency QPOs in the intermediate 
state moves to higher frequency and the low-frequency QPOs disappear when the state changes 
to HSS (see Belloni et al. 2006). 

These observations indicate that QPOs are associated with the cooling of the disk. By 
carrying out global three-dimensional (3D) magneto-hydrodynamic (MHD) simulations includ- 
ing radiative cooling, Machida et al. (2006) showed that when the accretion rate exceeds the 
limit for the onset of the cooling instability, the radiatively inefficient, optically thin disk tran- 
sitions into a magnetically supported, cool, intermediate state. Oda et al. (2007) constructed a 
steady model of such disks and showed that their luminosity can exceed O.lLEdd, where LEdd is 
the Eddington luminosity. When the transition to the cool disk takes place in the outer region, 
cool gas will be supplied to the inner region. 

Low frequency QPOs are sometimes attributed to the oscillation at the interface between 
the hot inner disk and the cool outer disk (e.g., Di Matteo & Psaltis 1999). The excitation 
mechanism of such an oscillation, however, was not clear. Another puzzle is the coexistence 
of low-frequency and high-frequency QPOs. Abramowicz & Kulzniak (2001) proposed that 
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high-frequency QPOs are generated by the resonance between radial and vertical oscillations 
of accretion disks. S. Kato (2001a,b) pointed out that non-axisymmetric g-mode oscillations 
can be trapped in a thin, relativistic disk and these oscillations are excited by the corotation 
resonance. Non-linear couplings of disk oscillations and disk warp were examined by S. Kato 
(2004). 

Y. Kato (2004b) reproduced high-frequency QPOs by 3D MHD simulations of optically 
thin, geometrically thick accretion flows. They showed that the QPOs appear around the 
epicyclic frequency n at 4rs and + /t, where VL is the Kepler frequency and showed that the 
amplitude of QPOs are damped after a while. 

The 1 //-noise-like fluctuations observed in BHCs have been reproduced by global 3D 
MHD simulations (e.g., Kawaguchi et al. 2000, Hawley 2001). Kawaguchi et al. (2000) showed 
that the PSD changes its slope around 10 Hz when they assumed a IOMq black hole. Machida 
& Matsumoto (2003) pointed out that magnetic reconnection in the innermost region of the 
disk can be the origin of intermittent X-ray flares known as X-ray shots (Negoro et al. 1995), 
which produce flat PSD at low-frequency (~ 1 Hz). 

In this paper, we report the results of global 3D MHD simulations which produced 
low-frequency QPOs and discuss their excitation mechanisms. 

In section 2, we describe basic equations and initial conditions. The results of simulations 
are given in section 3. In section 4, we concentrate on the oscillation excited in the inner torus. 
Section 5 is devoted to discussion and conclusion. 
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2. Numerical Methods 



2.1. Basic Equations 

We solved the following resistive MHD equations in a cylindrical coordinate system 



dp 



, V-(p«) = 0, (1) 



P 



dv ^ 



-VP-pV0 + ^^, (2) 



dB ^ , 47r 

— = Vx{vxB-—r^j) , (3) 

= ' 

where p, P, 0, -B, j = cV x B /An, r], T, and S" are the density, pressure, gravitational 
potential, velocity, magnetic field, current density, resistivity, temperature and specific entropy, 
respectively. The specific entropy is expressed as 5* = C^\n{P / p'^) , where is the specific heat 
capacity and 7 is the specific heat ratio. We included the Joule heating term but neglected 
the radiative cooling term in the energy equation. We assume the anomalous resistivity 77 = 
?7o[max(fd/fc — 1, 0)]^ (Yokoyama & Shibata 1994), where =j/p is the electron-ion drift 
speed and Vc is the threshold above which anomalous resistivity sets in. 

General relativistic effects are simulated using the pseudo-Newtonian potential = 
—GM/{r — Tg) (Paczyhski & Witta 1980), where G is the gravitational constant, M is the 
mass of the black hole, r = (gj^ + z'^Y^'^, and is the Schwarzschild radius. We neglect the 
self-gravity of the disk. 

2.2. Numerical Methods and Boundary Conditions 

We solved the resistive MHD equations using a modified Lax-Wendroff scheme (Rubin 
& Burstein 1967) with an artificial viscosity (Richtmyer & Morton 1967). 

The units of length and velocity are the Schwarzschild radius Vg and the light speed c, 
respectively. The unit time is to = rsC~^ = lO'^M/lOM© s. The unit temperature is given by 
To = mpC^^B^ = 1.1 X 10^3 K, where nip is the proton mass and is the Boltzmann constant. 

The number of grids is {Nt:^, N^p, Nz) = (250,64,384). The grid size is Aw = Az = 
0.1 for < w/vs < 10, and \z\/rs < 10. For model LT, we set the grid interval as follows; 
Awn = min(1.05Acan_i, Awmax) , Az-^ = min(1.05A2;n-i, Azmax), where Awma.x = lOAro and 
A^max = IOA2;. The outer boundaries ai vu = 132rs and ai z = ±70rs are free boundaries 
where waves can be transmitted. For model HT, we set Aw^^^ = 20Act7 and A^max = lOOA^. 
Therefore, the outer boundaries are located at w = 230rs and z = ±170rs, respectively. The 
grid size in the azimuthal direction is A(/9 = 27r/63. 

We included the full circle (0 < (/? < 27r) in the simulation region, and applied periodic 
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boundary conditions in the azimuthal direction. An absorbing boundary condition is imposed 
at r = Tin = 2rs by introducing a damping factor, 

D = 0.1 1.0-tanh ^ . (5) 

The physical quantities q = {p,v,B,P) inside r = ri^ are re-evaluated by 

g-w = g-D(g-go) , (6) 

which means that the deviation from initial values go is artificially reduced with damping rate 
D. Waves propagating inside r = are absorbed in the transition region (rin — 5 Aw <r < r^^). 

2.3. Initial Model 

The initial state of our simulation is an equilibrium torus threaded by a weak toroidal 
magnetic field. At the initial state, the torus is assumed to have a specific angular momentum, 
L oc w"-. 

The magnetic field distribution is determined according to Okada et al. (1989). By using 
the polytropic relation P = Kp"* at the initial state and by assuming 
^ SvrP ^ /t^7\ -2(7-1) 

where /3b is the initial plasma [3 at the initial pressure maximum of the torus {w,z) = {zuh,0), 
and is the azimuthal magnetic field. We integrated the equation of motion into a potential 
form, 

1 7 

*(cj7, 2;) = + — — H -v"^ + — —vl = ^h = constant , (8) 

7 — 1 2(7 — 1) 

where Vs = {'yP/pY^'^ is the sound speed, va = B^/ [ATipY/"^ is the Alfven speed, and \l/b = 
\l/(wb,0). At = tUb, the rotation speed of the torus L/w\, equals the Keplerian velocity. By 
using equation (8), we obtain the density distribution as 

[ max[vl^b-0-i.V(2t^^),O] Y'^'-'^ 

' 1.17/(7 -i)][i+ /3b- V2(.-)/^?-^)]/ ' 

where pb is the density at {w,z) = {zuh,0). Outside the torus, we assumed a hot, isothermal 
(T = Thaio) spherical halo. The density distribution of the halo is given by ph = Phaioexp[— (0 — 
0b)/(^B^haio)]5 where 0b is the gravitational potential at {zu,z) = {zub,0). 

In this paper, we report the results of simulations for two models. Model HT assumes 
a hot outer torus with the sound speed Cb = 0.029c at zu^ = SOrg. In model HT, the torus 
has a constant specific angular momentum (a = 0). Model LT is a cool disk model in which 
Cb = 0.01c at ti7b = 35rs and a = 0.43. In both models, we adopted /5b = 100, 7 = 5/3, L = 
{w]:,/2Y^'^zu^/ (wk — l)w", phalo = lO-^b, ^70 = 5 X 10"^, and Vc = 0.9c. Since we do not include 
radiative cooling, pb is arbitrary. We adopt pb = 1. 
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Fig. 1. (a) Time evolution of the mass accretion rate for model HT measured at tzj = 2.5rs. (b)Time 
evolution of the mass accretion rate at w = 2.5rs (black) and m = SOrg (gray), (c) Mass outflow rate 
measured at z = 40rs . 



Numerical Results 



3.1. A hot accretion disk : model HT 

Figure la shows the time evolution of mass accretion rate for model HT measured at 
w = 2.5rs. Figure lb shows the time evolution of mass accretion rate aX w = 2.5rs (black) and 
at G7 = 30rs (gray). The mass accretion rate M aX w = 2.5rs is computed by 

M = / pwv^dipdz . (10) 

J-20 JO 

The mass accretion rate at to = SOrg is measured by integrating the accretion rate in the 
equatorial region \z\ < Srg. The unit of the mass accretion rate is Mq = 27rp\^r^c. Mass accretion 
rate at to = 2.5rs becomes quasi-steady after about 10 rotational periods at the initial density 
maximum. The increase in mass accretion rate saturates when M ~ O.SMq. The mass accretion 
takes place due to the efficient angular momentum transport by Maxwell stress in an MRI-driven 
turbulent field. The equatorial mass accretion rate at to = SOrg shows time variation with a 
timescale of ~ 4000to- Figure Ic shows the time evolution of mass outfiow rate measured at 
z = 40rs computed by 

Mont = / pv^zudipdw . (11) 

J2 Jo 

The mass outfiow rate correlates with the mass accretion rate with a time lag of about 4000to- 
Figure 2a shows the correlation between the mass accretion rate at to = 2.5rs and the 
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Fig. 2. (a) Correlation between the mass accretion rate at tu = 2.5rs and the mass outflow rate at z = 40rs. 
(b) Self correlation of the mass outflow rate at z = 40rs. 

mass outflow rate measured at z = AOvg for model HT. The correlation is computed by 

F{6t) = f'° [{A{t) -A) -{3(1 + 6t) - B)]dt . (12) 

Here, A and B are average during the interval ts <t < t^- In figure 2, we adopted tg = ISOOOto 
and tg = 25000to- The correlation function has peaks at 6t = 3400to and 6t = 4800to- The mass 
outfiow rate at z = AOvs correlates with the mass accretion rate at zu = 2.5rs with delay of 
SOOOto — SOOOto- This timescale is comparable to the time scale of the propagation time of the 
outfiow whose average speed is ~ 0.01c. Figure 2b shows the self correlation function of the 
mass outfiow rate at 2; = 40rs in the time interval 21500 < t/to < 33500. In addition to the peak 
at St = 0, positive peaks appear around St = — 2400to and St = 2400to- This indicates that mass 
outfiow rate oscillates quasi-periodically with period ~ 2400to- 

Figure 3 shows the isosurface of the density and the vertical velocity. Blue surfaces 
and green surfaces depict the density isosurface (p = 0.2) and the isosurface of vertical velocity 
= 0.05c), respectively. Winds emerge intermittently from the inner region of the accretion 
disk. Toroidal magnetic fields are dominant in the wind. The intermittent ejection is driven by 
expansion of magnetic loops anchored to the accretion disk (Y. Kato et al. 2004a). 

Figure 4a shows the time evolution of magnetic energy integrated in 4 < ta/rs < 10, 
and \z\/rs < 1 normalized by the initial gas pressure at {zu,z) = {wh,0). Figure 4b shows 
the time evolution of 5^/87r averaged in 4 < cij/rs < 10, and \z\/rs < 1. Figures 4c and 4d 
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Fig. 3. Emergence of outflows from the inner region of the disk. Blue surfaces show the isosurface of the 
density {p = 0.2) for model HT at t = 29500io- Green surfaces show the isosurface of the vertical velocity 
(v^ = 0.05c). 
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Fig. 4. Time evolution of physical quantities for model HT averaged in 4 < vj/rs < 10, \z\/rs < 1, and 
< tp < 2tt. (a) Magnetic energy, (b) (B'^/Stt) / P\^, (c) time evolution of angular momentum transport 
rate, and (d) time evolution of the plasma /3. 
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Fig. 5. (a) Time development of the mass accretion rate AI ai w = 2.5rs for model LT. (b) Time evolution 
of the mass accretion rate M at nj = 2.5rs (black) and the equatorial mass accretion rate ai w = SOrg 
(gray). The interval 53000 < t/tQ < 70000 is enlarged, (c) Mass outflow rate at z = 40rs for model LT. 

show the time evolution of a = {B^B^/ A-k) / (P) and the ratio of the gas pressure to magnetic 
pressure /? = -PgasZ-Pmag, respectively. Although magnetic energy decreases, the plasma P is 
nearly constant (/3 ~ 8) because gas pressure also decreases. The fluctuations of the magnetic 
energy and radial magnetic field correlate with the mass accretion rate. 

3.2. A cool accretion disk: model LT 

Figure 5a shows the time evolution of the mass accretion rate for model LT measured 
at tz7 = 2.5rs. Figure 5b shows the mass accretion rate at w = 2.5rs (black) and the equatorial 
mass accretion rate at ccj = 30rs (gray) which enlarges the time range 53000 < t/to < 70000. 
Figure 5c shows the time evolution of the mass outflow rate measured at z = AOrg. In Figure 
5a, M ~ 0.04Mo is an order of magnitude smaller than that for model HT. This result indicates 
that the angular momentum transport rate strongly depends on the temperature of the gas 
supplied from the outer region. 

Figure 6a shows the correlation between the mass accretion rate at w = 2.5rs and the 
mass outflow rate at z = 40rs in the interval 52000 <t/tQ< 62000 for model LT. The mass outflow 
rate correlates with the mass accretion rate with time delay of 5000to — 6000to- Figure 6b shows 
the self correlation of the mass outflow rate at z = AOts in the interval 56650 < t/to < 68650. 
Positive peaks appear at 6t/to = —1500,0,2400,3500, and 4700. It indicates that the mass 
outflow rate oscillates with period lOOOto — 2000to- 
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Fig. 6. (a) Correlation between the mass accretion rate aX w = 2.5rs and mass outflow rate at z = 40rs. 
(b) Self correlation of the mass outflow rate at z = 40rs. 
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Fig. 7. Time evolution of (a) magnetic energy, (b) radial magnetic field, (c) angular momentum transport 
rate and (d) plasma /3 for model LT. 
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Fig. 8. The snapshots of the density distribution for model HT (left) and model LT (right), (a) and (b) 
show the density distribution averaged in the azimuthal direction, (c) and (d) show the density averaged 
in |z|/rs < 1. 

Figure 7 shows the time evolution of (a) magnetic energy, (b) (i?^/87r)/Pb, (c) angular 
momentum transport rate a, and (d) plasma (5 for model LT. The magnetic energy is averaged 
in the region 4 < w/rg < 10, \z/rs\ < 1, and < ip < 27t. 

Magnetic energy gradually decreases when t > 55000to- Figure 7c shows that a ~ 0.01 
in model LT. Since magnetic energy decreases, plasma f3 increases and stays around (3 ~ 10. 
Physical quantities shown in Figure 7 show short time scale oscillations and longer time scale 
(~ 4000to) time variations. The latter is due to the time variation of the mass accretion rate 
from the outer region. 

4. Formation of an Inner Torus and Its Oscillations 



Figure 8 shows the snap shots of the density distribution. Figures 8a and 8b show the 
density distribution in w — z plane for model HT and LT, respectively. The density is averaged 
in the azimuthal direction. Figures 8c and 8d show the density distribution in w — ip plane 
averaged in \z\/rs < 1. In model LT, an inner torus is created around w/r^ ~ 4 — 8. The 
inner torus is deformed into a crescent-like shape. The inner torus is formed because angular 
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Fig. 9. Radial distribution of physical quantities averaged in l^l/rs < 1 and in the azimuthal direction. 
Black curves show the distribution for model HT averaged in 23000 < i/io < 29000. Gray curves show the 
distribution for model LT averaged in 55000 < t/to < 61000. (a) Time-averaged specific angular momentum 
distribution. The dotted curve shows the Keplerian specific angular momentum using a Pseudo-Newtonian 
potential, (b) Radial distribution of the angular momentum transport rate defined by a= (B^B^/Att) / (P) . 
(c) Solid curves and dashed curves show the sound speed and Alfven speed, respectively, (d) Distribution 
of radial velocity. 

momentum transport becomes inefficient. Figures 7a and 7b indicate that Maxwell stress 
decreases due to the decrease in magnetic energy. Magnetic energy decreases partly because 
magnetic flux is swallowed into the black hole with accreting gas and partly because magnetic 
energy dissipates by magnetic reconnection. The deformation of the inner torus into a crescent 
shape takes place due to the growth of the Papaloizou-Pringle instability (Papaloizou & Pringle 
1984, Drury 1985). 

When the disk is hot, since the inner torus is not formed, disk gas accretes to the 
black hole through dense, spiral channels (Figure 8c). This result is consistent with the result 
obtained by simulations of the hot disk reported by Machida & Matsumoto (2003). 

Figure 9 shows the time-averaged radial distribution of physical quantities. The black 
curves show the results for model HT averaged in 23000 < t/to < 29000. The gray curve shows the 
result for model LT averaged in 55000 < t/to < 61000. Figure 9a displays the specific angular 
momentum distribution. The dotted curve shows the Keplerian specific angular momentum 
using pseudo-Newtonian potential. Since angular momentum continuously decreases in model 
HT, no inner torus is formed in this model. In model LT (dashed curve), since a is small, 
nearly constant angular momentum, inner torus is formed in oj/rg ~ 4 — 8. 
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Fig. 10. (a) and (b) w — z slice of distribution of temperature log(T/To) averaged in the azimuthal 
direction, (c) and (d) distribution of the vertical velocity averaged in the azimuthal direction, (a) and (c) 
are at t = 26950to for model HT, and (b) and (d) are at t = 62750to for model LT, 

Figure 9b shows the radial distribution of a = {B^B^pj An\) j (P) . Inside the inner torus, 
the angular momentum transport rate becomes very small. Solid and dashed curves in Figure 
9c show the sound speed and Alfven speed, respectively. Figure 9d shows the radial velocity. 
The accretion proceeds subsonically in > Srg. In model LT, accretion speed becomes very 
low in the inner torus. 

Figure 10 shows the distribution of azimuthally averaged temperature (10a, 10b), and 
vertical velocity (10c, lOd). Figures 10a and 10b show that accreting matter is significantly 
heated up. Figures 10c and lOd show that mass outflow emerges from the disk with an average 
speed of 0.01c. The outflow is anti-symmetric to the equatorial plane. The outflow is more 
powerful in model HT. By comparing Figures 8a and 10a, we can distinguish three regions; 
equatorial disk, hot funnel near the rotation axis, and outflows between the disk and the 
funnel. Accretion proceeds in the equatorial disk. A fraction of the accreting matter is ejected 
from the inner region of the disk. 

Figure 11 shows the distribution of the azimuthal magnetic field (color) and magnetic 
fields depicted from the poloidal components of magnetic fields averaged in the azimuthal 
direction (curves). The azimuthal magnetic field is antisymmetric to the equatorial plane at 
this stage (t = 26950to for model HT and t = 62750to for model LT). Magnetic fields are turbulent 
inside the disk but show more coherent structures in the interface between the disk and the 
halo where mass outflow takes place. Magnetic field lines are stretched along this interface. A 
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Fig. 11. Distribution of azimuthal magnetic field (color) and magnetic field lines (blue curves). The 
magnetic field lines are depicted from the poloidal components of magnetic fields averaged in the azimuthal 
direction, (a) model HT and (b) model LT . 

large-scale poloidal magnetic field is created in the funnel near the rotation axis. 

Figure 12 shows the snapshots of the density distribution on the vo — (f plane for model 
LT. Density is averaged in the same region as that in Figures 8c and 8d. The inner torus 
deforms its shape from a circle into a crescent, and from a crescent into a circle, repeatedly. 

The top panel of Figure 13 plots the time evolution of the ratio of the Joule heating 
rate to the magnetic energy averaged in 4 < w/r^ < 10, I^I/ts < 1, and < <y5 < 2ti. The ratio 
increases when magnetic energy is released (when (B'^/Stt) decreases or 77 increases). The 
ratio changes quasi-periodically with a time scale of ~ lOOOto- 

The bottom panel of Figure 13 shows the time evolution of the amplitude of non- 
axisymmetric m = 1 mode (m is the mode number in the azimuthal direction) of the density. 
The amplitude of the azimuthal mode is computed by Fourier decomposing the density contrast 



The running average of the amplitude is plotted using the amplitude during 250to for each 
point. The gray curve depicts the same curve as that in Figure 13. The amplitude of m = 1 
mode anti-correlates with rjj^ / {B"^ /Stt) . This indicates that the magnetic energy is released 
when the m = 1 mode disappears. The amplitude of the m=l mode also correlates with (-B^) 
(Figure 7b) and a (Figure 7c). They all show long timescale [t ~ 4000to) variation and short 
timescale {t ~ lOOOto) oscillations. The peaks of the short timescale oscillations in the Fourier 
amplitude of the m = 1 mode in Figure 13 coincide with those in Figure 7. 

Figure 14 shows the mass accretion rate measured at w = 4rs (black), magnetic energy 
(gray), and the Fourier amplitude of the m = 1 mode for the density (dashed) for model LT. 
When the magnetic energy increases inside the inner torus, the mass accretion rate increases 
because the angular momentum transport rate increases as magnetic energy (and magnetic 
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Fig. 12. Snapshots of the density distribution on ru — ip plane for model LT. 

stress) is accumulated in the torus. On the other hand, when the magnetic energy is released, 
the mass accretion rate decreases. The amplitude of the m = 1 mode (dashed curve in Figure 
14) correlates with the magnetic energy (gray). The magnetic fields are amplified when the 
amplitude of the m = 1 mode increases. On the other hand, the amplitude of the m = 1 mode 
decreases when the magnetic energy is released. When the Maxwell stress decreases due to the 
decrease in magnetic energy, the Papaloizou & Pringle instability grows again inside the torus. 
Therefore, the torus deforms itself into a crescent shape. The interval between the magnetic 
energy releases in the inner torus is ~ lOOOto- 

Figure 15a shows the correlation between the amplitude of the m = 1 mode and the Joule 
heating rate normalized by the magnetic energy for model HT. The correlation function has 
positive peaks at 6t/to = —2500,0 and 3000. The positive correlation ai 6t = indicates that 
magnetic energy is released when m = 1 spiral channel (see Figure 8c) develops. This result is 
consistent with that reported in Machida & Matsumoto (2003), where we showed that magnetic 
reconnection takes place in the spiral channel. Other peaks in the correlation function indicate 
that 771 = 1 mode develops quasi-periodically with interval 2500 — 3000to- 
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Fig. 13. (top) Time evolution of the ratio of the Joule heating rate to the magnetic energy for model 
LT. (bottom) Time evolution of the Fourier amplitude of the non-axisymmetric mode with the azimuthal 
mode number m = 1 computed from the density distribution (black). Gray curve depicts the same curve 
as that in the top panel. 

Figure 15b shows the correlation between the amphtude of the m = 1 mode and the Joule 
heating rate for model LT. In contrast to that in model HT, no positive peak appears at 5t = 0. 
Instead, negative peak appears at 5t/tQ = —1200 and positive peaks appear at St/to = 1200,2400. 
It indicates that m = 1 mode anti-correlates with the Joule heating rate and that the growth 
of the m = l mode precedes the release of the magnetic energy by St ~ 1200to- 

Figure 16 shows the spatial distribution of the Fourier amplitude z/P^ of time variabilities 
in mass accretion rate for model HT in 23000 < t/tQ < 32000, and for model LT in 55000 < 
t/tQ < 64000. In model HT, various peaks appear at various radius. On the other hand, in 
model LT, low frequency QPOs around lOHz appear in 5 < w/vs < 10, where the inner torus is 
formed. 

Figure 17 shows the Power Spectral Density (PSD) of the time variation of the mass 
accretion rate averaged in 3 < ta/vs < 8 and \z\/rs < 1. Black and gray curves show PSD {uPi, 
where is the Fourier power) for model LT and model HT, respectively. We adapted the data 
in the time range 23000 <t/to< 32000 in model HT and 55000 <t/to< 64000 in model LT. 

The PSD for model LT has a broad low-frequency peak around lOHz. This low-frequency 
peak corresponds to the oscillation involving amplification and release of magnetic energy in 
the inner torus. The PSD for model HT is fiat in u < lOHz, which means that oc u"^ and 
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Fig. 14. Correlation of the time evolution of the mass accretion rate and magnetic energy in model LT. 
(black) Time evolution of the mass accretion rate measured sX w = ^r^. (gray) Time evolution of the 
magnetic energy integrated in 4 < vo jr^ < 10, \z\/rs < 1, and < if < 2tt. (dashed) Fourier amplitude of 
the m = 1 mode in the density distribution. The curves are arbitrarily shifted in the vertical direction. 

changes its slope around u ~ 50Hz. Since oscillations are excited at various radius in model 
HT, PSD shows flat spectrum. The PSD in model LT has a slope steeper than that in model 
HT in 30Hz <u< lOOHz. 

In model LT in Figure 17, a small peak appears at z/ ~ 150Hz. Based on the simulation 
using 32 azimuthal grid points, Matsumoto & Machida (2007) pointed out that a high-frequency 
QPO appears when the disk shows low-frequency (~ 8Hz) sawtooth-like oscillations of magnetic 
energy. Here we confirmed that low-frequency QPOs and high-frequency QPOs coexist by a 
simulation including twice as many grid points in the azimuthal direction. 

5. Discussion 

In this paper, we studied the dependence of the structure and time variation of the 
black hole accretion flows on the gas temperature supplied from the outer region. When hot 
gas (T ~ 10^^ K) is supplied, the angular momentum is transported efficiently all the way to the 
black hole. The average ratio of the Maxwell stress to gas pressure, a = (B^B^/An) / (P) ~ 0.05, 
consistent with previous simulations by Hawley (2001) and Machida & Matsumoto (2003). Such 
disks correspond to the radiatively inefficient, optically thin disk in the low/hard state of black 
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Fig. 15. Correlation function between the amplitude of the m = 1 mode and the Joule heating rate 
normalized by the magnetic energy, (a) is model HT and (b) shows model LT. 
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Fig. 16. Radial distribution of the PSD of the time variation of mass accretion rate for model HT mea- 
sured in 23000 <t/to< 32000 (left) and for model LT measured in 55000 < t/to < 64000 (right). 
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Fig. 17. Power spectrum vPy, where is the Fourier power of the time variation of mass accretion rate 
averaged in 2.5 < vj/ts < 29 and \z\/rs < 1 for model LT (black) and for model HT (gray). 

hole candidates. Our simulation clearly showed the appearance of outflows with a maximum 
speed of ~ 0.05c from such hot disks. The mass outflow rate Mout is comparable to the mass 
accretion rate to the black hole. The magnetic fields are turbulent inside the disk but they 
show large-scale coherent structures near the rotation axis. Plasma flows out in the interface 
between the funnel and the accretion disk. The hot accretion flow shows time variations whose 
PSD is flat in < lOHz when the mass of the central black hole is M ~ IOMq. 

When cool plasma is supplied, an inner torus is formed around tJ7 ~ 4 — Srg. Such a 
torus is formed when the MRI-generated turbulent magnetic field is dissipated inside the disk. 
Since Maxwell stress decreases due to the decrease of magnetic energy, a decreases to a < 0.01. 
Thus a nearly constant-angular-momentum torus is formed in the innermost region {w < lOrg) 
of the accretion flow. We found that the inner torus deforms itself into a crescent shape. Such 
deformation takes place due to the growth of the non-axisymmetric instability in geometrically 
thick tori (Papaloizou & Pringle 1984). 

Such a non-axisymmetric structure enhances the growth of MRI. As magnetic energy 
increases, the angular momentum transport rate increases. Thus, the accretion rate increases. 
When the magnetic energy accumulated in the disk is released, the disk comes back to the 
weakly magnetized, axisymmetric torus. New cycle begins as magnetic energy is amplified 
by MRI. The period of the cycle is about lOOOto in model LT. It creates a low-frequency peak 
around 4 — 8Hz in PSD of the mass accretion rate. Low-frequency QPOs sometimes observed in 
low/hard state and hard intermediate state (HIMS) of black hole candidates can be reproduced 
by such magnetic cycles. We also showed that a high-frequency QPO is excited when low- 
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frequency QPO appears. 

Let us discuss why the magnetic cycle is excited when low-temperature gas is supplied. 
When hot gas accretes, since the scale height of the gas is large, large eddies are formed in the 
magnetically turbulent disk and create a large scale coherent magnetic field which transports 
angular momentum efficiently. On the other hand, when low-temperature gas is supplied, since 
the scale height and eddy size become smaller, angular momentum is transported only in the 
local region. Thus the angular momentum is transported more efficiently in high temperature 
disks. 

Low temperature plasmas can be supplied when cooling instability takes place in the 
outer disk. Machida et al. (2006) carried out a global 3D MHD simulation of black hole 
accretion flows by including radiative cooling. They showed that when the density of the outer 
disk is sufficiently high, cooling instability takes place. The outer disk shrinks vertically, and 
forms a magnetically supported, cool, optically thin disk. When such cool plasma accretes, the 
magnetic cycle may be excited in the inner region. 

The amplitude of the m = 1 non-axisymmetric mode correlates positively with the Joule 
heating rate in model HT but anti-correlates in model LT. In model HT, since the angular 
momentum transport rate is large enough, mass smoothly accretes along the spiral channel, in 
which the magnetic fields are stretched and form a current sheet. Magnetic energy is released 
in such current sheets. Therefore, the amplitude of the m = 1 mode correlates with the Joule 
heating rate. On the other hand, in model LT, non-axisymmetric m = 1 pattern disappears 
when the magnetic energy is released. 

Homan et al. (2005) pointed out that during the transition from the LHS to HSS, 
GX339— 4 shows a sub-transition from HIMS dominated by power-law X-ray radiation to soft 
intermediated state (SIMS) dominated by radiation from an optically thick disk. They also 
showed that the low-frequency QPOs appear when the X-ray spectrum stays in HIMS and 
SIMS and sometimes in LHS. The high-frequency QPOs appear in the HIMS and SIMS. When 
the high-frequency QPO is observed, a low-frequency QPO is always observed in the X-ray 
spectrum. This tendency is consistent with our simulation results. 

A number of theoretical models have been proposed for high-frequency QPOs. 
Abramowicz & Kulzniak (2001) proposed that high-frequency QPOs are formed by the res- 
onance of radial and vertical oscillations. S. Kato (2001a) studied the excitation of high- 
frequency QPOs by resonance with the disk warp. In our simulation, a high-frequency peak 
appears around u ~ 150Hz when the low-frequency oscillation is prominent. In our simulation, 
however, the low-frequency oscillation was weakened due to the heating of the inner torus. We 
expect that when the extra heating is extracted by radiative cooling, the magnetic cycle may 
continue and excite high-frequency QPOs. We would like to report the results of numerical 
simulations including radiative cooling in subsequent papers. 

In this work, we treated the relativistic effect by using the pseudo-Newtonian potential. 
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The accuracy of this approximation is worst near the black hole's horizon especially for spinning 
black holes. When we include the effects of black hole spin, the inner torus will be formed 
closer to the black hole. Thus, higher frequency QPOs will appear due to the oscillation of 
the inner torus. We should note that disk luminosity depends on the relativistic beaming and 
light bending. In subsequent papers, we would like to report PSDs obtained by relativistic 
ray-tracing of numerical results. 

We are grateful to M.A. Abramowicz, S. Kato, W. Kulzniak, and M. Burusa for dis- 
cussion. Numerical computations were carried out on VPP5000 at Center for Computational 
Astrophysics, CfCA of NAOJ (P.I. MM). This work is supported in part by Japan Society for 
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